% By Martin M. Andreasen, Auguest 2009
% This function computes the conditional moments up to k_period into the
% future for the variables selected by "y_select" in the g-function. 
% Note that y_select must have the form of 1 x dimy, where dimy is the
% number of variables in the g-functions we want to compute conditional
% moments for. The steady state value of y_select should be stored in y_ss
%
% In terms of the notation we solve: p(x,sig) = E_t[r(x_t+1,sig)]
% The computations are derived by using the Perturbation On Perturbation (POP) method.
% This function is for a log-transformation
%
function [px,pxx,pss,pxxx,pssx] = CondMoments_log(gx,gxx,gss,gxxx,gssx,hx,hxx,hss,hxxx,hssx,eta,y_select,k_period);

% Allocating memory
dimy    = size(y_select,2);
nx      = size(hx,1);
ne      = size(eta,2);
px      = zeros(dimy,k_period,nx);
pxx     = zeros(dimy,k_period,nx,nx);
pss     = zeros(dimy,k_period);
pxxx    = zeros(dimy,k_period,nx,nx,nx);    
pssx    = zeros(dimy,k_period,nx);


for i=1:dimy
    rx   = reshape(gx(y_select(1,i),:),1,nx);
    rxx  = reshape(gxx(y_select(1,i),:,:),nx,nx);
    rss  = gss(y_select(1,i),1);
    rxxx = reshape(gxxx(y_select(1,i),:,:,:),nx,nx,nx);
    rssx = reshape(gssx(y_select(1,i),:),1,nx);
    
   
    for j=1:k_period
       % ************** The first order effects *****************
       px(i,j,:) = rx(1,:)*hx;
       
       
       %*************** The second order effects ***************
       tmp = zeros(nx,nx);
       for gama1=1:nx
           tmp = tmp + rx(1,gama1)*squeeze(hxx(gama1,:,:));
       end
       pxx(i,j,:,:) = hx'*rxx*hx + tmp;
      
       pss(i,j) = rx(1,:)*eta*eta'*rx(1,:)' + trace(eta'*rxx*eta) + rx(1,:)*hss + rss;
       
       
       % ***************** The third order effects *****************
       for alfa1=1:nx
           for alfa2=1:nx
                tmp = zeros(1,nx);
                for gama3=1:nx
                    tmp = tmp + hx(:,alfa1)'*reshape(rxxx(:,:,gama3),nx,nx)*hx(:,alfa2)*hx(gama3,:);
                end
                pxxx(i,j,alfa1,alfa2,:) = tmp  ...
                                        + hx(:,alfa1)'*rxx*reshape(hxx(:,alfa2,:),nx,nx) ...
                                        + hx(:,alfa2)'*rxx*reshape(hxx(:,alfa1,:),nx,nx)...
                                        + reshape(hxx(:,alfa1,alfa2),1,nx)*rxx*hx(:,:)...
                                        + rx(1,:)*reshape(hxxx(:,alfa1,alfa2,:),nx,nx);
           end
       end       
       
       pssx(i,j,:) = 2*rx(1,:)*eta*eta'*rxx*hx+hss'*rxx*hx+rx(1,:)*hssx+rssx*hx;
       for gama3=1:nx
           pssx(i,j,:) = reshape(pssx(i,j,:),1,nx) + trace(eta'*squeeze(rxxx(:,:,gama3))*eta)*hx(gama3,:);
       end
       
       % Updating rx
       rx   = reshape(px(i,j,:),1,nx);
       rxx  = reshape(pxx(i,j,:,:),nx,nx);
       rss  = pss(i,j);
       rxxx = reshape(pxxx(i,j,:,:,:),nx,nx,nx);
       rssx = reshape(pssx(i,j,:),1,nx);       
    end
end

% We reexpress the output if dimy = 1
if dimy == 1
    px      = squeeze(px(dimy,1:k_period,1:nx));
    pxx     = squeeze(pxx(dimy,1:k_period,1:nx,1:nx));
    pss     = reshape(pss(dimy,1:k_period),k_period,1);
    pxxx    = squeeze(pxxx(dimy,1:k_period,1:nx,1:nx,1:nx));    
    pssx    = squeeze(pssx(dimy,1:k_period,1:nx));
end